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ABSTRACT 

Fermi has provided the largest sample of 7-ray selected blazars to date. In 
this work we use a complete sample of FSRQs detected during the first year of 
operation to determine the luminosity function (LF) and its evolution with cosmic 
time. The number density of FSRQs grows dramatically up to redshift ~0. 5-2.0 
and declines thereafter. The redshift of the peak in the density is luminosity 
dependent, with more luminous sources peaking at earlier times; thus the LF of 
7-ray FSRQs follows a luminosity-dependent density evolution similarly to that 
of radio-quiet AGN. Also using data from the Swift Burst Alert Telescope we 
derive the average spectral energy distribution of FSRQs in the 10 keV-100 GeV 
band and show that there is no correlation of the peak 7-ray luminosity with 
7-ray peak frequency. The coupling of the SED and LF allows us to predict that 
the contribution of FSRQs to the Fermi isotropic 7-ray background is 9.3ti % 
(±3% systematic uncertainty) in the 0.1-100 GeV band. Finally we determine 
the LF of unbeamed FSRQs, finding that FSRQs have an average Lorentz factor 
of 7 = II.7I2' 3 ,, that most are seen within 5° of the jet axis, and that they 
represent only ~0.1 % of the parent population. 
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Introduction 



The detection of luminous quasars at redshift >6 (e.g. iFan et al.l 120031 ; IWillott et al. 



20101 ) provides evidence of super-massive black hole (SMBHs) forma tion in the first 1 Gyr 



of cosmic time. There are appreciable challenges to forming (s ee e.g. I Wyithe &: Loeb 



Volonteri fc Reesl 20051 iBegelman et alJliooS iMaver etllhoioh and fueling fseelKauffmann Sz Haehnelt 



2003 



2000l : IWvithe fc Loebl 20031 ; Croton et al.ll2Q06h these objects at such early ti mes, although it 



is widely believed that strong accretion can be in itiated by major mergers (see lKauffmann fc Haehnelt 
2OO0t IWvithe fc Loebl 120031 : ICroton et al.ll2006h . 



Blazars represent an extreme manifestation of such AGN activity, with radiation along 
the Earth line-of-sight dominated by a relativistic jet. It is, as yet, unclear how such jet 
activity conn ects with the more iso t ropic ally emitted bulk accretion luminosity. For example 



according to iBlandford fc Znajekl ( 119771 ). the energy stored in a black hole's spin can be 
extracted in the form of a relativistic jet. Thus blazar evolution may be connected with 
the cosmic evolution of the spin states of massive black holes. Radio- loud (jet dominated) 



blazars have been seen at redshifts as high as z=5.5 (IRomanil 120061 ). and it is plausible 



Escala et al. 


2004; 


Dotti et al. 


2007) 



Thus the study of radio-loud (RL) AGN, blazars, with strong relativistically beamed 
jets can provide a method to study jet activity, BH spin, and major merger events. This 
can be done by determining the luminosity function of blazars (LF, essentially the number 
of blazars per comoving volume element within a certain luminosity range) and its evolution 
with redshift. The Fermi Gamma-ray Space Telescope provides one of the largest data 
sets with which to study the properties of blazars. Thanks to its sensitivity and uniform 
coverage of the sky, Fe rmi has detected hundreds of blazars from low redshifts out to z=3.1 
(lAbdo A.. etaDboilh . 



The LF of blazars also allows us to evaluate their contribution to the diffuse back- 



grounds and to determine their relationship with the parent popul ation (lAiello et al. 



Inoue 



201 



Wall et al. 



2001 



2005), soft X-ray ( Giommi fc Padoyani 19941 Rector et 



Caccianiga et al.ll2002 



Beckmann et al 



2003; 



]. Blazars have been extensively studied at radio ( Dunlop fc Peacock _ 199ol: 



Padovani et al 



2000 



2008b; 



Wolter & Celotti 



20071 ) and GeV energies 



( IHartman et al.lll999l ). It seems clear that flat spectrum radio quasars (FSRQs) evolve posi- 
tively (i.e. there were more blazars in t he pastjDunlop &; Peacock 119901 ) up to a redshift cut - 



off which depends on luminosity (e.g. iPadovani et al. 



2007 



Wall 2008 



Aiello et al.ll2009bh . 



In this respect FSRQs evolve similarly to the population of X-r ay selected, radio-quiet, AGNs 
( jUeda et al.l 120031 : lHasinger et al.ll2005t lLa Franca et al.l 120051 ) . On the other hand, the evo- 
lution of the other major class of Fermi-detected AGN, BL Lac objects, and their relation 



-3 - 



Padovani et al.l 120071 ) or even negative evolution (e.g. 



no evolution i 


Caccianiea et al. 


2002; 


Rector et al. 


2000: 


Beckmann et al. 



20031 ). Samples with larger redshift completeness fractions are needed to study the LF of 



these claims. 



In this work we report on the LF of FSRQs dete cted by Fermi in it s first year of 
operation. There have been attempts in the past (e.g. Chiang et aL 1995b to character - 



ize the evolution of 7-ray AGN, starting from the EGRET sample (IHartman et al.lll999l ). 
One challenge was the small sample size and redshift incompleteness of the EGRET set. 
Often it was assumed that the 7-ray detected blazars had a LF follow ing that of another 



„ — - - .! ~~ 

(1996). Narumoto & Totani 


(2006). Inoue & Totani (200< 


)). Stecker & Venters (2010) fol- 


low this approach. Alternatively, a LF may be estimated from the 7-ray sample directly 


(Chiane & Mukheriee 1998: 


Miicke & Pohl 2000: Dermer 


2007: Bhattacharva et al. 2009). 



although only a small (~60) blazar sample including both BL Lac objects and FSRQs was 
available (from EGRET data) with acceptable incompleteness. 

We report here on a detailed LF measured from a sample of 186 7-ray selected FSRQs 
detected by Fermi. The work is organized as follows. Sections [2] and [3] describe the properties 
of the sample used and the method employed to determine the LF of blazars. The luminosity 
function of FSRQs is derived in Section [U In Section [5] the spectral energy distributions 
(SEDs) of Fermis FSRQs are analyzed in detail, testing for possible correlations of the 
peak energy with the peak l uminosity. The con tribution of FSRQs to the isotropic gamma- 
ray backgroun cfl(IGRB, see lAbdo et al.l l2010dl ) is determined and discussed in Section |6j 
Throughout this paper, we assume a standard concordance cosmology (H =71 km s _1 Mpc -1 , 
ft M =l-ft A =0.27). 



2. The Sample 



The First Fermi LAT Catalog (1FGL, lAbdo et al.l l2010al ) reports on more than 1400 
sources detected by Fermi- LAT during its first year of operation. The first LAT AGN catalog 
(1LAC, lAbdo et al.l hoioh associates ~700 of the high-latitude 1FGL sources (|b| > 10°) 
with AGN of various types, most of which are blazars. T he sample u s ed for this analysis 
consists of sources detected by the pipeline developed by lAbdo et al.l (|2010eJ) with a test 



1 The isotropic gamma-ray background refers to the isotropic co mponent of the Fer mi sky ( Abdo et al 



2010dJ) and as such might include components generated locally (e.g. iKeshet et al.ll2004l ) and components of 
truly extragalactic origin. 



-4- 



statistic (TS) significance greater than 50 and with |b| >15°. For these sample cuts we have 
produced a set of Monte Carlo simulations that can be used to determine and account for 
the selection effects. This sample contains 483 objects, 186 of which are classified as FSRQs. 
The faintest identified FSRQ has a 100 MeV - 100 GeV band flux of F 100 « 10" 8 ph cm" 2 
s _1 . To limit the incompleteness (i.e. the fraction of sources without an association) we 
apply this as a flux limit to the sample of 483 objects, resulting in a full sample of 433 
sources of which 29 (i.e. ~7%) do not have associated counterparts. 

The composition of this sample is reported in Table [TJ The 186 FSRQs detected by 
Fermi with TS> 50, |b| > 15°, and Fi o > 10 -8 ph cm -2 s _1 constitute the sample that will 
be used in this analysis. 



Method 



A classical approach to derive the LF is based on the 1 /Vmax method of lSchmidtl (119681 ) 
applied to redshift bins. However, this method is known to introduce bias if there is sig- 
nificant evolution within the bins. Moreover, given our relatively small sample size and 
the large volume and luminosity range spanned, binning would result in a loss of informa- 
tio n. Thus we de c ided to apply the maximu m-likelihood ( ML) al gorithm first introduced 
by iMarshall et al.l ( 119831 ) and used recently by lAjello et al.l ( I2009bl ) for the study of blazars 
detected by Swift. The aim of this analysis is to determine the space density of FSRQs as 
a function of rest-frame 0.1-100 GeV luminosity (L 7 ), redshift (z) and photon index (r), by 
fitting to the functional form 



d 3 N 



dL^dzdT 



d 2 N dN dV 

x x — 

dL^dV dT dz 



_ r , dN dV 



(i) 



where $(L 7 , z) is the luminosity function , and dV/d z is the co- moving volume element per 



unit redshift and unit solid angle (see e.g. lHogpl999l ). The function dN/dT is the (intrinsic) 
photon index distribution and is assumed to be independent of z. It is modeled as a Gaussian: 



dN 
~dY 



(r-M) z 



(2) 



where fi and a are respectively the mean and the dispersion of the Gaussian distribution. 

The best-fit LF is found by comparing, through a maximum-likelihood estimator, the 
number of expected objects (for a given model LF) to the observed number while accounting 
for selection effects in the survey. In this method, the space of luminosity, redshift, and 
photon index is divided into small intervals of size dL 1 dzdT. In each element, the expected 
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Table 1. Composition of the |b| >15°, TS>50, Fi o > 10 8 ph cm 2 s 1 sample used in 

this analysis. 



CLASS 


# objects 


Total 


433 


FSRQs 


186 


BL Lacs 


157 


Pulsars 


28 


Other a 


16 


Radio Associations b 


17 


Unassociated sources 


29 



a Includes starburst galaxies, LIN- 
ERS, narrow line Seyfert 1 objects 
and Seyfert galaxy candidates. 

b Fermi sources with a radio coun- 
terpart, but no optical type or red- 
shift measurement. 
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number of blazars with luminosity L 7 , redshift z and photon index V is: 

A(L 7 , z, T)dL~ j dzdT = $(L 7 , z)Q(L 7 , z, r)— — -— ^-dL 7 dzdT (3) 

where Q(L 7 ,z, T) is the sky coverage and represents the probability of detecting in this 
survey a blazar with luminosity L 7 , redshift z and pho ton index T. This probability was 
derived for the sample used here by lAbdo et al.l (l2010el ) and the reader is referred to that 
aforementioned paper for more details. With sufficiently fine sampling of the L 1 — z — V 
space the infinitesimal element will either contain or 1 FSRQs. In this regime one has a 
likelihood function based on joint Poisson probabilities: 

L = Y[ A(L 7)i , Zi, T i )dL 7 dzdTe- x ^ i ' Zi ' Ti ^ dzdr x JJ e -A(i 7J -,*j,r J )dt r d»dr ^ 

» 3 

This is the combined probability of observing one blazar in each bin of (L 7) j, z iy I\) populated 
by one Fermi FSRQ and zero FSRQs for all other (L 7i j, Zj, Tj). Transforming to the standard 
expression S = —2 In L and dropping terms which are not model dependent, we obtain: 

The limits of integration of Eq.[5j unless otherwise stated, are: L ljrnin =10 u erg s _1 , L 7 max =10 52 erg 
s _1 , z min =0, z max =6, T min =1.8 and T max =3.0. The best-fit parameters are determined by 
minimizingG S and the associated 1 a error are computed by varying the parameter of inter- 
est, while the others are allowed to float, until an increment of AS—1 is achieyed. T his gives 
an estimate of the 68% confidence region for the parameter of interest (lAvnil 119761 ) . While 
computationally intensive, Eq. [5] has the advantage that each source has its appropriate 
individual detection efficiency and k-correction treated independently. 

In order to test whether the best-fit LF provides a good description of the data we 
compare the observed redshift, luminosity, index and source count distributions against the 
prediction of the LF. The first three distributions can be obtained from the LF as: 

dN 

dz 

" -■- min 

dN <- r 



dL 7 
dN 

~d~r 



\{L 1 ,T,z)dL 1 dT (6) 

—"v min 
i a x 

\(L 7 ,r,z)dzdT (7) 

71 

I I \(L 7 ,r,z)dL 7 dz (8) 

** L^v .m.i.n. ** Zm.im 



L. 



2 The MINUIT minimization package, embedded in ROOT (root.cern.ch), has been used for this purpose. 



-7- 



where the extremes of integrations are the same as in Eq. [5j The source count distribution 
can be derived as : 



N(> S) 



.dNdV, vl 
$(L 7 ^) _ . dTdzdL^ 



dT dz 



(9) 



where L^(z, S) is the luminosity of a source at redshift z having a flux of S. 



To disp lay the LF w e rely on the ''N obs /N method devised bv lLa Franca fc Cristiani 



(11997) andlMivaii et al.l (20011 ) and employed in several recent works (e.g. 



La Franca et al 



20051 ; lHasinger et al.ll2005l ). Once a best-fit function for the LF has been found, it is possible 



to determine the value of the observed LF in a given bin of luminosity and redshift: 

N obs 



^*(-^ / 7,i) 



J\jrndl 



(10) 



where and z% are the luminosity and redshift of the i th bin, $ md/ (L 7i j, Z{) is the best- 
fit LF model and N° bs and j\T mdZ a re the observed and th e predicted number of FSRQs in 
that bin. These two techniques (the iMarshall et al.l (jl983h ML method and the "N ofcs /N mrf '" 



estim ator) provide a minimally biased estimate of the luminosity function, (cf. iMiyaji et al. 
200lh . 



4. Results 

4.1. Pure Luminosity Evolution and the Evidence for a Redshift Peak 

The space density of radio-quiet AGNs is known to be maximal at inte rmediate redshift . 



The epoch of 1 


;his 'r 


Hasineer et al. 


2005) 



cosmic time and a fall-off in fueling activity as the rate of major mergers decreases at late 
times. To test whether such behavior is also typical of the LAT FSRQ population, we 
perform a fit to the data using a pure luminosity evolution (PLE) model of the form: 



$(L 7)2 ) = $(L 7 /e(z)) 



where 



and 



<$>{L,/e{z = 0)) 



dN 



A 



dL 7 ln(10)L 7 



e(z) = 1 + zYe 



72 



-i -1 



(11) 

(12) 
(13) 
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In this model the evolution is entirely in luminosity: i.e. the FSRQ were more luminous 
in the past if positive evolution (k > 0) is found (the opposite is true otherwise). It is also 
straightforward to demonstrate that the luminosity evolution (i.e. Eq. ITB"]) of FSRQs peaks 
at z c =— 1 — k£. The best-fit parameters are reported in Table [2J The evolution of the FSRQ 
class is found to be positive and fast (k = 5.70 ± 1.02). The redshift peak is z c = 1.62 ±0.03. 
Moreover, the subsequent rate of decrease of the luminosity after the peak is well constrained 
(£ = —0.46 ± 0.01). However, as shown in Fig. [TJ while this model provides a good fit to 
the observed redshift and luminosity distributions, it is a very poor representation of the 
measured logAMogS'. 




Fig. 1. — Redshift (left), luminosity (middle) and source count (right) distribution of LAT 
FSRQs. The dashed line is the best-fit PLE model discussed in the text. 

We next test whether the redshift peak depends on luminosity, splitting the data set 
at L 7 = 3.2 x 10 47 erg s _1 . A fit is then performed to each sub-sample to determine the 
position of the redshift peak (if any), keeping the parameters of Eq. [T2l fixed. The results of 
the fits to the low- and high-luminosity data sets are reported in Table [2j From Eq. [13] and 
the values of k and £ it is apparent that there is a significant shift in the redshift peak, with 
the low- and high- luminosity samples peaking at ~1.15 and ~1 .77, respectively. 

Table 2. Best-fit parameters of the Pure Luminosity Evolution LF. Parameters without 
an error estimate were kept fixed during the fit. 



Sample 


# Objects 


A a 


71 


L* 


72 


k 


t 


M 


a 


ALL 


186 


5.59(±0.41)xl0 3 


0.29±0.53 


0.026±0.066 


1.25±0.32 


5.70±1.02 


-0.46±0.01 


2.45±0.13 


0.18±0.01 


Low L 


89 


15.4(±0.2)xl0 3 


0.29 


0.026 


1.25 


4.30±2.39 


-0.50±0.04 


2.47±0.04 


0.21±0.03 


High L 


97 


22.6(±2.0)xl0 3 


0.29 


0.026 


1.25 


3.47±1.73 


-0.79±0.04 


2.46±0.02 


0.20±0.02 



a In unit of 10 13 Mpc 3 erg 1 s. 
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4.2. The Luminosity Dependent Density Evolution and the Redshift Peak 



Since a simple PLE LF model provides an inadequate fit to the Fermi data and since 
there is some evidence for the evolution of the redshift peak with luminosity, we now fit 
the Fermi FSRQ set to a luminosity-dependent density evolution model (LDDE). Here the 
evolution is primarily in density with a luminosity-dependent redshift peak. The LDDE 
model is parametrized 



where 



and 



e(z, L„ 



7> ■ 



$(L 7 ) x e(z, Lj) 



l + z 



1 + z c {L^ 



pi 



l + z 



1 + z c (L 7 ) 



P 2 



z c (L 7 ) = z* c ■ (LJ10 



48\a 



(14) 



(15) 



(16) 



$(L 7 ) is the sa me double power law used in Eq. [12j This parametrization is similar 



to that proposed by lUeda et al.l (120031 ). but is continuous around the redshift peak z c (L 7J 
This has obvious advantages for fitting algorithms that rely on the derivatives of the fitting 
function to find the minimum. Here z c (L^) corresponds to the (luminosity-dependent) red- 
shift where the evolution changes sign (positive to negative), with z* being the redshift peak 
for a FSRQ with a luminosity of 10 48 erg s _1 . 

The LDDE model provides a good fit to the LAT data and is able to reproduce the 
observed distribution in Fig. |2j The log-likelihood ratio test shows that the improvement 
over the best PLE model is significant, with a chance probability of ~ 10 -6 . Results are 
reported in Table El 

In Fig. [3] we subdivide the sample into four redshift bins with comparable number of 
sources to illustrate how the LF changes. The evolution, visible as a shifting of the peak 
and a change of the shape of the LF between different bins of redshift, is clearly visible. 
This evolution takes place mostly below redshift ~1.1 where the space density of our least 
luminous objects (i.e. L 7 rs 10 46 erg s _1 ) increases by ~ lOx. Above this redshift the 
variation is less marked, but one notices that: 



• the space density of logL=47 objects decreases from redshift 1 to redshift 1.5 while 
that of logL=48 FSRQs still increases (lower left panel, green versus red line). This 
indicates that the space density of logL=47 FSRQs peaks at a redshift ~ 1.1. 

• similar behavior holds for logL=48 FSRQs in the highest redshift bin so that their 
maximum space density should occur well below z=3. 
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The best-fit parameters confirm that the redshift of maximum space density increases 
with increasing luminosity (with the power-law index of the redshift-peak evolution a = 
0.21 ± 0.03, see Eq. [To]) . This redshift evolution can be clearly seen in Fig. HJ which shows 
the change in space density for different luminosity classes. 




Fig. 2. — Redshift (upper left), luminosity (upper right), photon index (lower left), and 
source count (lower right) distributions of LAT FSRQs. The dashed line is the best-fit 
LDDE model convolved with the selection effects of Fermi. Notice the greatly improved 
source count distribution over the predictions of the PLE model of Figure [TJ 



Table 3. Best-fit parameters of the Luminosity Dependent Density Evolution LF 



Sample 


# Objects 


A a 


71 


L* 


72 




a 


Pi 


P2 


M 


a 


ALL 

ALL b 

ALL C 


186 
208 
186 


3.06(±0.23)xl0 4 
2.82(±0.19)xl0 4 
8.72(±0.63)xl0 3 


Q.21±0.12 
0.26±0.12 
0.38±0.16 


0.84±0.49 
0.87±0.53 
0.89±0.70 


1.58±0.27 
1.60±0.27 
1.60±0.30 


1.47±0.16 
1.42±0.15 
1.38±0.18 


0.21±0.03 
0.20±0.03 
0.18±0.03 


7.35±1.74 
8.21±1.78 
7.71±1.84 


-6.51±1.97 
-5.66±1.73 
-4.44±1.78 


2.44±0.01 
2.42±0.01 


0.18±0.01 
0.19±0.01 



a In unit of 10 -13 Mpc -3 erg^ 1 s. 

b 22 unassociated sources were included in this sample by drawing random redshifts from the observed redshift distribution of FSRQs. 
c Derived using the detection efficiency for curved reported in Fig, 1161 
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Fig. 3. — LF of the Fermi FSRQs in different bins of redshift, reconstructed using the 
N bs/N mc ii method. The lines represent the best-fit LDDE model of § 14.21 To highlight the 
evolution, the LF from the next lower redshift bin is over-plotted (dashed lines). 

4.3. Analysis of Uncertainties 

One of the main uncertainties of our analysis is due to the incompleteness of the FSRQ 
sample. In Table HJ there are 17 sources with associated radio counterparts lacking optical 
type and redshift measurements. A fraction of these may be FSRQs. In addition, there 
are 29 sources without any statistically associated radio counterpart. The lack of radio flux 
means that these cannot be FSRQ similar to those in the Fermi sample, unless position 
errors have prevented radio counterpart associations. Thus even if a few of these sources are 
mis-localized, the maximum possible incompleteness of our FSRQ sample is on the order of 
20/(186+20), i.e. -10%. 



The standard way to account for this incompleteness is to correct upwards the normal- 



13 




Fig. 4. — Growth and evolution of different luminosity classes of FSRQs. Note that the space 
density of the most luminous FSRQs peaks earlier in the history of the Universe while the 
bulk of the population (i.e. the low luminosity objects) are more abundant at later times. 
The range of measured distribution is determined by requiring at least one source within the 
volume (lower left) and sensitivity limitations of Fermi (upper right). 



ization of the LF as to reflect the likely real number of FSRQs (associated or not) in our 
sample. Considering extra information about these sources, this likely incompleteness is even 
smaller than the 10% above. First, we find that only 50% of all the radio-loud identified 
sources in our sample are FSRQs. This suggests that only 8 of the 17 radio identified sources 
are FSRQs, with the remainder being BL Lac objects and lower luminosity AGN. A similar 
argument can be made based on the 7-ray spectral index. The median index for FSRQs is 
T = 2.44 while only 9% of our BL Lac objects have such soft spectra. There are 4 such 
soft sources in our set of 17 radio sources. Conservatively assuming that all are FSRQs, we 
infer that twice this number, namely 8, FSRQ are in the radio-detected s ample. Of the 29 



source s without radio counterparts, 4 are Blazar 'ANTI-Associations' (see lAbdo et al.ll2009 



2010j . for details), where we can definitively exclude any flat-spectrum radio source bright 
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enough for a Fermi-type blazar. While most of these 29 sources lack the very deep radio 
observations to make an ANTI-Association, all but 5 have been classified as pulsar candi- 
dates, based on gamma-ray spectral curvature and lack of variability We thus suspect that 
virtually all of these sources are other classes of gamma-ray emitters, eg. pulsars yet to be 
discovered, starburst galaxies, etc. Our conclusion is that the likely incompleteness is only 
8/(186+8) =4%. Conservatively adding a few nominally radio-quiet sources from erroneous 
LAT localizations may allow 5% incompleteness. This correction has been applied in § 14.11 
and §H2J 

Moreover, we show that our results are robust even in the unlikely case that half of 
the unassociated sources are FSRQs (thus neglecting that the unassociated FSRQs occupy a 
definite part of the parameter space). In this scenario, the number of likely yet unidentified 
FSRQs is ~22 (i.e. 8 + 29/2). We assign a random redshift, drawn from a smoothed version 
of the FSRQ redshift distribution to 22 out of the 46 unassociated objects (i.e. giving 10 % 
maximum incompleteness for the FSRQ sample). These are then used with the associated 
FSRQs to derive the luminosity function. One example is reported in Table [31 It is apparent 
that including the unassociated sources with random redshifts does not modify the shape 
and the parameters of the LF. Indeed, all the parameters are well within the statistical 
uncertainties of the parameters of the LF derived using only the associated data set. This 
test shows that the systematic uncertainty introduced by the incompleteness is smaller than 
the statistical uncertainty. In Appendix [A] we present a more detailed discussion of other 
sources of uncertainty. 



4.4. Comparison with Previous Results 

4-4- 1- The Local Luminosity Function 



The local LF is the luminosity function at redshift zero. For an evolving population, the 
local LF is obtained by de-evolving the luminosities (or the den sities) ac c ordin g: to the best- 
fit model. This is generally done using the 1/Vmax method of ISchmidtl (119681 ). as reported 
for example by lDella Ceca et all (120081 ). However, since the best representation of the LF is 
the LDDE model, the maximum volume has to be weighted by the density evolution relative 
to the luminosity of the source. In this case, the maximum allowed volume for a given source 
is defined as: 

e(z,Li) dV 



Vmax 



£l(Li,z) 



&\Z'm,i'n.i Li) (Lz 



dz 



(17) 



where Lj is the source luminosity, Q(Li, z) is the sky coverage, z max is the redshift above which 
the source drops out of the survey, and e(z, Li) is the evolution term of Eq. [15] normalized 
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(through e(z min , Li)) at the redshift z min to which the LF is to be de-evolved. Th e LF de- 
evolved at z min (z min =0 in this case) is built using the standard 1/Vmax method (jSchmidt 



19681 ). 



In order to gauge the uncertainties that the different methods might introduce in the 
determination of the local LF we consider also an alternate method. We perform a Monte 
Carlo simulation, drawing 1000 series of parameters from the covariance matrix of the best 
fit LDDE model described in § 14.21 Using the covariance matrix ensures that parameters 
are drawn correctly, taking into account their correlations. The re-sampled parameters are 
used to compute the ±1 a error of the LF at redshift zero. This is reported in Fig. |5j There 
is very good agreement with the local LFs using this method and the 1/Vmax approach. 
The gray band in Fig. shows the true statistical uncertainty on the space density that the 
1/Vmax method (applied using the best-fit parameters) is not able to capture. 

We find a local LF described by a power law with index of 1.6-1.7, for Fioo < 10 47 erg 
s _1 , steepening at higher luminosity. Thus the local LF can be parametrized as a double 
power law: 



dN 
dL^ 



A 



71 



+ 



72 



where A = (3.99±0.30) x 10" 11 , L* = 0.22±0.30, 7i=1.68±0.17, 7 2 =3.15±0.63 and both L 7 
and L* are in units of 10 48 erg s _1 . Other models (e.g. a Schecter function, a simple power 
law etc.) do not generally provide as good a fit to the data. The values of the low-luminosity 
index ji and the high-luminosity i ndex 79 are in g ood a greement with that found here of 
1.63±0.16 and 2.3±0.3 reported bv lPadovani et all J2OO7I ) for the DRBXS survey of FSRQs. 



Values very similar t o those found here were also reported for a radio FSRQ sample by 
Dunlop fc Peacock ( 1990 ). who fincE^i = 1.83 and 72 = 2.96. The Fer mi LF low-luminosity 
i ndex (i.e. 71) is flatter than that determined using EGRET blazars by lNarumoto fc Totani 
(120061 ) as is app arent in Fig. \5[ How ever, a re-analysis of the same data sets employing the 
blazar sequence ( Fossati et al. 1998h to mo del the blazar SEDs found a lo w-luminosity index 
in the 1.8-2.1 range (llnoue fc Totanill2009l ). Also, in a more recent work, llnoue et al.l (120101 ) 
modified their SED models to be able to reproduce TeV data of known blazars. Their LF 
at redshift zero (see Fig. [5]) is found to be in relatively good agreement with that found here 
for the Fermi sample. 



3 Their definition of local luminosity function and Ecu 1181 differ by a 1/L, (or 1/P in their paper) term. 
Thus we added 1.0 to their exponents. 
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Fig. 5. — Local (z=0) LF of the Fermi FSRQs as derived from the best-fit LDDE model 
in § I4.2I (solid line). The gray band represents the ±1 a uncertainty computed as described 
in the text. The long - and short-dashed l i nes sh ow t he LF model s base d on the EGRET 
blazars derived by by iNarumoto &: Totanil (120061 ) and llnoue et al.l (I2010f)respective ly. The 
dashed-dotted line shows the prediction from the model of IStecker fc Venters! (120101 ). 



4-4- 2. The Luminosity Function at Redshift 1 

Fig. [6] shows the luminosity of FSRQ at redshift 1 compared to predictions from recent 
models. It is apparent that the evolution of the Fermi LF is stronger than predicted by 
any of these models. The increase in space density from redshift to 1 for a source with 
a luminosity of 10 48 erg s _1 is almost a factor ~ 150. This dramatic increase is not seen in 
the evolution of radio-quiet AGN (e.g. lUeda et al.l l2Q03t lHasinger et al.l 120051 ) whose space 
density increases by a factor 25 -50 between redshift and 1. The increase of a factor ~60 
seen in FSRQs detected in Radio iDunlop &: Peacock! ( I1990l ). is still slower than that of Fermi 
blazars. This explains why the predictions based on luminosity functions derived at other 
wavelengths (see Fig. [6]) underpredict the density of high-luminosity Fermi FSRQs at redshift 
of 1. 
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Fig. 6. — LF of the Fermi FSRQs at redshift 1.0 as derived from the best-fit LDDE model in 
§ 14.21 (solid line). The gray band represents the ±1 a uncertainty computed with the method 
described in the text. The long dashed and sh ort dashed lines show the LF models based 



on the EGRET blazars derived respectively by lNarumoto fc Totani] (120061) and llnoue et al. 



(2010). The dashed-dotted line shows the prediction from the model of IStecker fc Venters 
J2010h . 



5. The Spectral Energy Distributions of FSRQs 

We may use the 0.1-300 GeV LAT spectra of our uniform bright Fermi FSRQ sample 
along with the 15-200 keV spectra measured by the Swift-BAT to characterize the high 
energy Inverse Compton (IC) sector of the blazar SED. In turn this allows us to describe 
the average SED properties of the energetic blazars and to estimate their contribution to the 
cosmic high-energy backgrounds. 
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5.1. Data Analysis 



For this analysis, we further restrict the sample to Fioo > 7 x 10 _8 ph cm -2 s" 1 (cor- 
responding to an energy flux of 3.4xl0 _11 erg s" 1 cm" 2 for a power law with an index of 
2.45), since for brighter sources Fer mi has a negligible bias in the detected spectral indices 
and FSRQs are selected uniformly (jAbdo et al.ll2010el ). For fainter sources hard-spectrum 
objects (principally BL Lac objects) domi nate the sample. There are 103 bright FSRQs 
detected by Fermi that meet these criteria ( jAbdo et al.ll2010il ). 

We analyze two years of Fermi data using version V9r21 of the science toolfl The 
data are filtered, removing time periods in which the instrument was not in sky-survey 
mode and photons whose zenith angle is larger than 100 degrees. We consider only photons 
collected within 15 degrees of the source position with 100 MeV<E< 300 GeV. We employ 
the P7SOURCE_V6 instrumental response function (IRF) and perform binned likelihood 
analysis using the gtlike tool. First, a likelihood fit using a power-law model for all the 
sources in the region of interest (ROI) is performed on the entire energy band (100 Me V- 
300 GeV). The parameters (i.e. flux and photon index) of all the sources within 3° of the 
target FSRQ, along with the normalization of the diffuse mod el, are left f ree. More distant 
sources have parameters frozen at the 2FGL measured values (jAbdoll201ll ). We next choose 
30 logarithmically spaced energy bins and perform a binned likelihood in each, deriving the 
flux of the target FSRQ in each energy bin. During this exercise the flux of the FSRQ is 
allowed to vary while the photon index is fixed at the best-fit value found for the whole 
band. All the neighboring sources had parameters fixed at the best-fit values, although 
the diffuse emission normalization was allowed to vary. This analysis provides a 30-bin 
100 MeV-300 GeV energy spectrum for all 103 sources in the bright FSRQ sample. 

Swift-BAT is a coded- mask telescope that has conducted a seve ral year survey in the 
15-200 keV hard X-ray sky tjGehrels et al.ll2004l ; iBarthelmv et al.ll2005l ). With this deep expo- 
sure, BAT reaches a sensitivity of ^10 ~ n erg cm" 2 s" 1 on most of the high-latitude sky (e.g. 
Tueller et al]l2007t kjello et al.ll2008al lcl ICusumano et alJboioh. Blazars rep resent 15-20% 
of the extragalactic source population detected by BAT (lAjello et al.ll2009bl ). Since Fermi- 
LAT and Swift-BAT have comparable sensitivity in their respective bands^ and since the two 
bands cover the bulk of the IC component, a joint study allows an accurate characterization 
of the IC spectrum and the contribution to the background. 



4 http://fermi. gsfc.nasa.gov/ssc/data/analysis/software/ 



5 A Fermi FSRQ with a photon flux of 3 x 10" 8 ph cm" 2 s" 1 in the 100 MeV-100 GeV band and a power-law 
spectrum with an index of 2.4 has a energy flux ofl.5x!0 -11 erg cm" 2 s" 1 . 
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We use ~6 years of BAT data to extract a 15-200 keV spectrum for all the FSRQs in 
the Fermi s a mple. Spectral extraction is perf ormed accor d ing to the r ecipes presented b y 
Aiello et all J2008ch and discussed in detail bv lAiello et all ((2009a]) and kiello et all feoioh . 
Both BAT and LAT spectra are multiplied by 4nDi(z) 2 (with Dl(z) the luminosity distance 
at redshift z) to transform the flux into a luminosity and shifted by (1+z) to transform into 
source rest-frame SEDs. 



For each FSRQ, we fit the BAT and LAT data with an empirical model of the following 



form: 



9 dN , n9 



E 



E 



1BAT / 



1LAT' 



. e -^WK . A-kD l {z) 2 (19) 

where Jbat and Jlat are the power- law indices in the BAT and the LAT bands and and 
E c are the break and the cut-off energy, respectively. The e~^/ E l Ec term allows us to model 
the curvature that is clearly visible in a few of the Fermi spectra. The fit is performed only 
for E< 20 GeV to avoid possible steeping due to the absorption of 7-ray photons by the 
extra-galactic background light (EBL; e.g. Stecker et al. 2006 ; Franceschini et allboOo ). 



Two sample spectra are shown in Fig. [7J It is apparent that for the brightest FSRQs, 
BAT and LAT are efficient in constraining the shape of the IC emission. Even when the 
BAT signal is not significant, the upper limit from BAT still provides useful constraints on 
the low energy curvature of the SED. In a number of bright sources (see e.g. Fig. [7J the 
highest-energy datapoint in BAT at >120keV is seen to deviate from the baseline fit. This 
deviation is at present not significant (i.e. the reduced x 2 of the baseline fit is already ~1.0), 
but certainly suggestive of a second component. Observations with INTEGRAL extending 
to energies >200 keV might ascertain the nature of this feature. 

Several caveats necessarily apply to our analysis. First, the BAT and LAT observations 
are not strictly simultaneous. LAT spectra are accumulated over 2 years while the BAT 
data span 6 years. In principle one could restrict the BAT data to the period spanned by 
the Fermi observations. In practice this would seriously limit the BAT sensitivity, weaken- 
ing constraints on most of the spectra. Second, it is possible that BAT and LAT are not 
sampling exactly the same emission component. In particular, BAT might be domi nated by 
IC emission produced by the synchrotron self Compton (SSC, iMaraschi et al.lll992l) compo- 



nent while the LAT may be more sensitive to External Compton (EC. lDermer fc Schlickeiser 
19931 ) emission. Ultimately detailed SED modeling with strictly simultaneous data would be 
needed to eliminate these concerns, and such work is well beyond the scope of this paper. 
Bearing these caveats in mind, our bright sample is nearly free of selection effects other than 
the hard flux-limit threshold applied to the Fermi data. This allows an detailed study of the 
average properties of the high-energy SED of FSRQs. 
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Two BAT-LAT spectra of famous blazars fitted with the empirical model described 



5.2. Correlation of Peak Luminosity and the Energy of the Peak 

We can compare the Inverse Compton rest-frame peak luminosity and peak energy 
from the SED fits to the FSRQ in our sample. As shown in Figure El there is no apparent 
correlation between these quantities; indeed the Kendall test gives a value r = 0.09 indicating 
no significant correlation. Since generally neither Fermi nor BAT directly sample the high- 
energy peak, we also fit the spectra using a third degree polynomial function instead of the 
model in Eq. UHl The results are shown in the right panel of Fig. [8] and confirm the previous 
findings. 



This is in contrast to the correlation found (but see also iNieppola et al.ll2008l) between 



the luminosity and the energy of the synchrotron peak of blazars (e.g. iGhisellini et al.lll998 
Fossati et al.lll998l ). This might imply that the jet parameters (e.g. Doppler factor, lumi- 
nosity of the target photon field, etc.) do not depend on blazar GeV luminosity or redshift. 
This may be understood if the IC peak is largely controlled by EC emission for these sources. 



5.3. Average SEDs 



It is useful to estimate the average SED of FSRQs, particularly for estimating the 
contribution of FSRQs to the extragalactic gamma-ray background. First we define the 
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Fig. 8. — Peak luminosity versus the energy of the peak for the complete sample of FSRQs 
discussed in § [5j The left plot shows the value derived using a double power law with expo- 
nential cut-off while the right panel show parameters derived using a third degree polynomial 
function. 

bolometric luminosity as the luminosity in the 10 keV - 300 GeV banc|§ and divide the sources 
into four bins of bolometric luminosity with approximately the same number of objects in 
each bin. In these luminosity bins we compute the average of the logarithm of the spectral 
luminosity at each energy. Associated errors on this average spectrum are computed using 
the Jackknife technique. In this framework we neglect uncertainties due to different level 
of the energy density of the extragalactic background light which would affect mostly the 
high-energy part of the SED (i.e. >20GeV). 

Fig. [9] shows the average rest-frame SED for the FSRQ sample in the four luminosity 
bins. This plot confirms the lack of a systematic correlation of the peak luminosity and 
energy. Indeed, all the averaged SEDs show a peak in the 10-100 MeV band and their shape 
does not change much with luminosity. 

To transform luminosities between observed and rest-frame we need the k-correction, 
along with its redshift variation, shown in Fig. [TUJ In practice, there is little difference 
bet ween the k-correction for the average SED computed here (even applying EBL absorption, 
e.g. lFranceschini et al.ll2008l ) and one computed for a power law (i.e. (l+z) r ~ 2 ) with a photon 
index of 2.4. Only at large redshifts do the two k-corrections start to differ; this difference is 
only ~5% at a redshift of 4. We find that using a power law index of 2.37 and taking into 



The best fit is extrapolated from 20 GeV to 300 GeV. 
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account EBL absorption allows us to reproduce correctly the k-correction up to redshift ~6. 
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Fig. 9. — Average rest-frame spectral energy distributions for four representative FSRQ 
luminosity classes (left panel), see § 15.31 In each SED, the band represents the la uncer- 
tainty on the average. This does not reflect the uncertainty connected to the level of the 
extragalactic background light. 



6. The Contribution to the Isotropic Gamma-Ray Background 



The nature of the diffuse gamma-ray background at GeV energies remains one of the 
most interesting problems in astroph ysics. The presence of a n isotropic component was first 
determined by the O SO-3 satellite (IKraushaar et al.l 119721) an d confirmed by SAS-2 and 
EGRET (respectively iFichtel et al.lll975t ISreekumar et al.lll998f ). This isotropic component 
is normally referred to as the isotropic gamma-ray background (IGRB). Fermi recently pro- 
vided a refined measurement of this isotropic component showing that it can be adequately 
described as a single power law with an index of 2.4 in the 200 MeV - 100 GeV energy range 
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Fig. 10. — Ratio of source rest-frame luminosity to observed luminosity (i.e. k-correction) 
as a function of redshift for the average SED shape reported in § 15.31 and for two generic 
power laws. 



(jAbdo et al.ll2010dh . 



Blazars, representing the most numerous identified populations by EGRET and Fermi 



dictions range from 20-30 % ( 


Chiang & Mukheriee 


1998; 


Miicke k Pohl 


2000; 


Narumoto & Totani 


2006 


Dermer 


2007; 


Inoue & Totanilhoom to 100 % 


Stecker & Salamon 


1996|; 


Stecker & Venters 



20101 ). Analysis of the source-count distribution showed that for F 100 > 10 ph cm s 
the contribution of unresolved blazars to the IGRB is ~16 % in the 100 MeV - 100 GeV band 
( jAbdo et al.ll2010ef ). Since the source counts distribution show a strong break at a flux of 
.Fioo ~ 6 x 10 _8 ph cm -2 s _1 , it was concluded that extrapolating the source counts to zero 
flux would produce ~23 % of the IGRB. 



Now, with a measured LF we can more robustly evaluate the emission arising from 
faint FSRQs. In addition, the FSRQ SED shape study of the previous section also allows 
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improvement over the simple power-law type spectra assumed by lAbdo et al.l (l2010el ). 



The contribution of 'unresolved' FSRQs to the IGRB can be estimated as: 



Fegb 



, dV 
dz—— 

dz 



dT 



dL 7 F 7 (L 7 , z) 



d 3 N 



dL y dzdT 



(20) 



where the limit of integration are the same as those of Eq. Eland F 1 (L 11 z) is the flux of 
a source with rest-frame lumino sity L 7 at redshift z. Since we are interested in the diffuse 
flux not yet resolved by Fermi (lAbdo et al.ll2010df ) the term (l-fl(r, -F 7 )/fi ma:r ) takes into 
acco unt the photon ind ex and source flux dependence of the LAT source detection threshold 



sec 



Abdo et al.l l2010d . for more details). The limits of integration of Eq. [20] are the same 



as those of Eq. El 

In Eq. [20] setting Q(T, F 7 ) /Vt max =Q allows us to compute the total 7-ray emission arising 
from the FSRQ class. The result is 3.13±^x lO" 6 ph cm" 2 s" 1 sr" 1 in the 100 MeV - 100 GeV 
band. This value should be compared with the tot al sky intensity of 1 .44xl0~ 5 ph cm -2 s _1 
sr -1 , which includes IGRB plus detected sources (lAbdo et al.ll2010dl ). Thus FSRQs make 



2l.7ti b 7 % of this total intensity. 



If one considers the contribution only from the FSRQs that Fermi has not detected then 
this becomes 9.66j£{g x 10" 7 ph cm 2 s 1 sr 1 (with a maximum systematic uncertainty of 

Al . T his represents 9.31} q% of the IGRB intensity in the 
2010dl ). From above it is also clear that Fermi has already 



sr 



see 



3x10 "ph cm 2 s 
0.1-100 GeV band ( lAbdo et al. 



resolved more than 50 % of the total flux arising from the FSRQ class. Fig. [IT] shows this 
contribution. The possible presence of external Compton components in the SEDs of FSRQs 
makes the estimate between 200 keV and 100 MeV uncertain (see § 15. 3p . Future observations 
with both Fermi above 20 MeV and INTEGRAL above 200 keV and physical modeling of 
blazar spectra might substantially reduce this uncertainty. 

Even the (disfavored) PLE model cannot accommodate a much larger contribution of 
FSRQs to the IGRB. Indeed, in this case the contribution of unresolved FSRQs would be 
1.37xl0~ 6 ph cm" 2 s" 1 sr" 1 (or ~13% of the IGRB intensity). 



7. Beaming: The Intrinsic Luminosity Function and the Parent Population 

The luminosities L defined in this work are apparent isotropic luminosities. Since the jet 
material is moving at relativistic speed (7 >1), the observed, Doppler boosted, luminosities 
are related to the intrinsic values by: 



L = 5 P 3? 



(21) 
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Fig. 11. — Contribution of unresolved (left) and all (right) FSRQs to the diffuse extragalactic 
background (blue line) as determined by integrating the luminosity function coupled to the 
SED model derived in § 15.31 The hatched band around the best-fit prediction shows the 1 a 
statistical uncertainty while the gray band represents the systematic uncertainty. 

where Jf is the intrinsic (unbeamed) luminosity and 5 is the kinematic Doppler factor 

5 = (j - vV - 1 cos 9) 1 (22) 

where 7 = (1 — /3 2 ) _1//2 is the Lorentz factor and /3 = v/c is the velocity of the emitting 
plasma. Assuming that the sources have a Lorentz factor 7 in the 71 < 7 < 72 range 
then the minimum Doppler factor is 5 m i n = 7^ (when #=90°) and the maximum is S max = 
72 + \Z~f2 ~ 1 (when 6 = 0°). We adopt a value of p = 4 that applies to the case of continuous 
jet emission which seems appropriate for the study presented here since long-term average 
luminosities are used. The case p = 4 applies also to spherical blobs if the observed emission 
is dominated by the SSC component , while a val ue of p=5-6 should be adopted if the 
emission is due to external Compton (lDermerlll995l ). However, as shown later these values 
imply extremely small isotropic rest-frame luminosities. 



Beaming is known to alter the shape of the intrinsic luminosity function. lUrry fc Shafer 



( 11984J ) provide an analytic solution to the case where the intri nsic luminosity function is a 
single power law and the jets have a single Lorentz factor. In lUrry fc Padovanil ( 1l99ll ) the 
intrinsic luminosity function may be a double power law and a distribution of Lorentz factor 
is considered. 



In this Section we will determine the intrinsic luminosity function of the Fermi FSRQs 
and their Lorentz and Doppler factor distributions. In what follows we adopt the formalism 
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and symbols of iListerl (120031 ) and ICara fc Lister! (120081 ). 

We begin by defining the intrinsic luminosity function as: 



(23) 



valid in the S£\ < jSf < j£?2 ra nge. The pro bability of observing a beamed luminosity L given 
a Doppler factor 5 is (see also IListerl 120031 ): 



P(L,6) = P S (6)*$(J?) 



7l 



(24) 



where is the probability density for the Doppler 5. Assuming a random distribution 

for the jet angles (i.e. Pg = sin 6), this results in 



From here it follows that 



P,{l)Pe{e) 



Ps(S) 



dO 



dS 



dj 



72 P 7 (7) 



(25) 



(26) 



where -P 7 (7) is the probability density for 7 and the lower limit of integ ration f(5) depends 
on the Doppler factor value and is reported in Eq. A6 in IListerl (120031 ) . Integrating over 5 
yields the observed luminosity function of the Doppler beamed FSRQs: 



$(L) = fciL 



—B 



S 2 (L) 



Ps(S)6 



p(B-l) 



5 1 (L) 



where, as in ICara fc Lister! ( 120081 ). the limits of integration are 

6x(L) = min{5 ma:E ,max (S min , (L/^ 2 ) llp )} 
5 2 {L) = max{5 min ,min (5 max , (L/^fi) 1/p )} 



(27) 



(28) 
(29) 



In this way, by fitting Eq. [27] to the Fermi Doppler boosted LF, it is possible to determine 
the parameters of the intrinsic luminosity function and of the Lorentz-factor distribution. 

We assume that the probability density distribution for 7 is a power law of the form 

P 7 ( 7 ) = (30) 

where C is a normalization constant and the function is valid for 71 < 7 < 72. Here we 
assu me 71 = 5 and 72 = 40 as this is the range of Lorentz factors observed for 7 -loud blazars 
(e.g. lLahteenmaki fe Valtaojall2003t iLister et al.ll2009at ISavolainen et al.ll20101 ). While the 



largest intrinsic luminosity can be set free, the lowest one depends on the value of p 
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chosen: i.e. from Eq. [21] the beaming factor defines the intrinsic luminosity corresponding to 
the apparent isotropic luminosity we observe. For p=4 and p=5, 5£\ has to be set 10 40 erg 
s _1 and 10 38 erg s -1 respectively. We set S£i = lO A J?i, but this choice has hardly any impact 
on the results. 

The free parameters of the problem are: the normalization (k\) and the slope (B) of 
the intrinsic LF and the slope k of the Lorentz factor distribution. We have fitted Eq. [27] to 
the Fermi LF de-evolved at redshift zero derived in §4.4.11 Fig. [T2] shows how the best-fit 
beaming model reproduces the local LF of FSRQs measured by Fermi. From the best-fit 
we derive, for the p = 4 case, an intrinsic LF slope of B — 3.04 ± 0.08 and an index of the 
Lorentz-factor distribution of k\ = —2.03 ± 0.70. The fit values are summarized in Tabled] 
The Lorentz-factor distribution implies an average Lorentz factor of de tected Fermi blazars o f 



7 = 11. 7lg 2> i n reasonable agreement with measured values (see e.g. iGhisellini et al.ll2009l ). 



The index of the Lorentz-factor distr ibution is in agreement with k\ ~-1.5 found for the 



CJ-F survey ( iLister fc Marscherlll997l ). The parameters for the p =5 case are very similar 
to those of the p =4 case, but the reduced \ 2 is slightly worse (see Tabled]). Nevertheless, 
the predictions of the two models are in agreement and we find again that the average bulk 
Lorentz factor is 7 = 10.2^ 2 8 . As noted already the extreme Doppler boosting (<5 5 ) requires 
the intrinsic luminosities to be small: i.e. 10 38 erg s _1 < ££ <10 42 erg s -1 . 

From the ratio of the integrals of the beamed and intrinsic LF we derive that the 
Fermi FSRQs represent only 0.1% of the parent population. The average space density 
of LAT FSRQs (derived from the LF § 14. 2p is 1.4 Gpc -3 , implying that the average space 
density of the parent population is ~1500Gpc~ 3 . Our model also allows us to determine 
the distribution of jet angles with respect to our line of sight. This is found to peak at 
~1.0 degrees (Fig. [13]). While FSRQs can still be detected at large (i.e. >10 degrees) off- 
axis angles for reasonably low 7 factors (~5-7), most FSRQs detected by Fermi are seen at 
angles less than 5 degrees from the jet axis. Owing to their larger space density (see Fig. [T2"]) 
misaligned jets produce a non-negligible diffuse emission. From our model it is found that 
the ratio between the diffuse emission contribution of misaligned jets and that of blazars 
(at redshift zero) is ~30%. This has obvious consequences for the generation of the IGRB. 
In fact nearly all of the flux produced by radio galaxies is unresolved, with only 2 steep- 
spectr um radio quasars, and 2 FR II and 7 FR I radio galaxies detected with the Fermi 



LAT ( lAbdo et al.l 1201 Obi ). All the results reported above apply to both the p = 4 and p = 5 



models. 
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Fig. 12. — Fermts LF de-evolved at redshift zero and the best-fit beaming model (for p=4, 
see text) described in § [71 The red continuous line shows the predicted space density of 
misaligned jets. 



8. Discussion and Conclusions 



In this paper we examine the properties of 7-ray selected FSRQs using data from the 
Fermi-LAT instrument. Our work relies on a nearly complete, flux-limited sample of 186 
FSRQs detected by Fermi at high significance and large Galactic latitude during the first 
year of operations. This analysis explores several of the properties of FSRQs; here we discuss 
and summarize our findings. 



8.1. Beamed Luminosity Function 



The redshift-zero LF of Fermi FSRQs is well described by a double power-law model, 
typical for the LF of AGN (both of the radio-quiet and radio- loud). At mid-to-high lumi- 
nosities th ere is good agreement between the Fermi LF and that determined using EGRET 
data (e.g. iNarumoto fc Totanil l2006t llnoue et al.ll20T0l ). At luminosities < 10 46 erg s -1 the 
FSRQ LF appears to be slightly flatter than in previous studies. 
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Fig. 13. — Distribution of viewing angles with respect to the jet axis for Fermi FSRQs. 



The space density of LAT-detected FSRQs increases dramatically with redshift, growing 
by 50-100 x by z—1.5. Describing the evolution of the LF as simple luminosity evolution 
(PLE model), there are strong indications that the evolution must cut off for z> 1.6. After 
this redshift, the space density of blazars starts to decrease quickly. 

A simple PLE does not fully explain the Fermi data. In particular the source count 
distribution is not well modeled. Since there is evidence that low- and high-luminosity 
sources have different redshift peaks, we consider a more sophisticated model where the 
evolution is primarily in density, but objects with different luminosity are allowed to have 
different redshift peaks. This so-called LDDE m odel explains well the evolutionary beh avior 



of (radio-quiet) AGN s elected in the X-ray band (jUeda et al.ll2003l ; lHasinger et al.ll2005f ) and 



was also suggested by iNarumoto fc Totanil (120061 ) to describe the LF of EGRET blazars 



■cp 

Stecker & Ventersll2010f ) are not in agreement with the LF of Fermi FSRQs at redshift unity. 
This is due to the fact that the Fermi FSRQs is found to evolve more quickly than the LFs 
of X-ray-selected AGN or radio-selected FSRQs. Indeed, the space density of Fermi FSRQs 
increases by a factor ~150 between redshift and 1 while the density increase is at most a 



Narumoto & Totani 


2006; 


Inoue et al. 


2010 
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factor ~60 for the models discussed above. 



The LDDE model implies that sources with a luminosity of 10 erg s 



erg s 

and 10 48 erg s _1 reach their maximum space density at a redshift of ~0.6, ~0.9, and ~1.5 
respectively. It is clear, then, that the most luminous objects, while lower in numbers, appear 
before the bulk of the (low-luminosity) population. This is the first time that this is seen in 7- 
rays. This "anti-hierarchica l" behavior, where the largest structure s come first is common to 
all classes of AGN (see e.g. ICowie et al.lll999t lHasinger et al.ll2005l . and references therein), 
and is often referred to as "cosmological downsizing". The disappearance of quasar-like 
objects at late times mi g ht indicate that accretio n efficiency ev olves as a function o f cosmic 
time (e.g. lMerlonill2004l ). |Pi Matteo et al.l ( 120051 ) (but see also ISanders et al.lll988l ) propose 
that the merging of two massive galaxies leads to, in addition to strong star formation 
activity, a burst of inflow feeding gas to the SMBH and initiates a 'quasar-like' phase. 
Eventually the energy released by the AGN in the form of powerful winds expels the gas, 
quenches star formation and starves the AGN. This pi cture, coupled with the fact that major 
merg ers become increasingly rare at low redshift (e.g. iFakhouri et al.ll2010l ; iKulkarni fc Loeb 
20 111 ) may explain why quasars are rare in the local Universe. 



Fig. [TH shows the energy density injected in the Universe (e.g. the luminosity density) 
by FSRQs as function of redshift. This shows a broad peak between a redshift of 1 and 2. 
A sim ilar behavior is shown by the cosmic star-formation history (e.g. iHopkins &: Beacom 



20061 ) which peaks around redshift 1-2. This represents a strong link between the host 
and the nucleus. A noteworthy fact is the mild evidence for a fast decline in the space 
density of FSRQs after the redshift peak (see parameter P2 in Table [3]). The decline seems 
to be as dramatic as the increase in space density leading up to the redshift peak. For 



comparison, X-ray sele cted samples of AGN show a much mi 



the reds hift peak (e.g. lUeda et al.l 120031 ; lHasinger et al 



2005 



der decline 



Aird et al. 



p2 



-1.5) after 
2010|). However, 



recently Silverman et al. ( 20081 ) (but see also Schmidt et al. 1995 ) reported evidence for a 



similarly dramatic decrease in the space density of AGN. 

One factor contributin g to this phenomenon is the difficulty for Fermi to detect soft 
sources (lAbdo et al.l l2010d ). At redshift >3 the SED peak should move well below the 
current LAT energy band, making it difficult to probe a population of extremely soft sources. 
Increasing the effective area at or below 100 MeV may help uncover such a population. 
Because the rising part of the IC peak is in the hard (>10keV) X-ray band, high-redshift 



objects are more easily selected in this band (see e.g. the Swift/BAT results in lAjello et al. 



( ]2009bl )). In this case another strategy would be to build a bolometric luminosity function 
that uses both the 7-ray and the X-ray selected samples. 
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Fig. 14. — Luminosity density as a function of redshift produced by the Fermi FSRQs. The 
gray band represents the 1 a statistical uncertainty around the best-fit LF model. 



8.2. The Intrinsic Luminosity Function 

Doppler boosting allows Fermi to detect many blazars when their jet emission is within a 
few degrees from the line of sight. As shown first by lUrry fc Shaferl (119841 ). Doppler boosting 
is known to alter the shape of the luminosity function. In this paper, we adopted a formalism 
that allowed us to recover the intrinsic de-beamed LF and to determine the distribution of 
bulk Lorentz factors for the Fermi FSRQs. 



The intrinsic LF is compatible with a single steep power law with an index of 3.04±0.08 
in the range of intrinsic luminosities 10 40 erg s _1 < Jz? < 10 44 erg s _1 . The break seen in 
the beamed LF at redshift zero is thus produced by Doppler boosting. The data cannot be 
explain by a single, averaged, Lorentz factor, but require a distribution of Lorentz factors. 
This distribution is found to be compatible with a power law with an index of -2.03±0.70 
in the 5<7<40 range. This yields the result that the average FSRQ bulk L orentz factor 
is T = ll.TtfJ, i n good agreement with several studies ( iGhisellini et al .112 0091 ). Our model 
is able to predict the distribution of viewing angles with respect to the jet axes of Fermi 
FSRQs. It is found, see Fig. [T3J that on average FSRQs are seen within an average angle of 
~2.3° from the jet and that most are seen within 5°-6°. A few Fermi FSRQ detections may 
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Fig. 15. — Number density of LAT-detected FSRQs as a function of redshift. The gray band 
represents the 1 a statistical uncertainty around the best-fit LF model. 



be up to 15° off-axis (if these have low Doppler factors). Fermi-detected FSRQs represent 
only ~0.1 % of the parent population for randomly pointed jets. 

Monitoring observations with the Very Long Baseline Array (VLBA) established that 
LAT detected bla zars have, on average, significantly faster ap parent jet speeds than non-LAT 
detected blazars (ILister et al.ll2009bl ; ISavolainen et al.l 120101 ) . Their distribution of Lorentz 
factors is in good agreement with the results of our analysis. Moreover, they report the 
distribution of viewing angles with respect to the jet axis for FSRQs detected by LAT. From 
their study the average viewing angle is 2.9° ± 0.3° and all the FSRQs in their sample have 
9 < 5°. There is thus substantial agreement within the VLBA monitoring observations and 
the results of our beaming model applied to 7-ray only data. 

The space density of FR-II radio galaxies (i.e. the putative parent popula tion of FSRQs) 
i s repo rted to be ~1580 Gpc" 3 f at 15 GHz) and ~2200Gpc" 3 (at 1.4 GHz) bv lCara fc Lister 



(120081 ) and iGendre et al.l (120101 ) respectively. From our study we derive a space density of 
FSRQ parents of ~1500Gpc -3 in substantial agreement with the numbers above. 



Future work may test whether the intrinsic properties of blazars (i.e. Lorentz factor, 
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luminosity etc.) evolve with redshift. This will likely require larger complete samples and 
improved modeling of selections effects. 



8.3. Spectral Energy Distribution 



Blazars SEDs are characterized by the typical "two hump" spectrum where the low- 
energy peak is produced by electrons radiating via synchrotron and the high energy peak 



scenario; 


Maraschi et al. 


1992) 


Dermer & Schlickeisei 


1993). 



In this work we have combined quasi-simultaneous Swift/BAT and Fermi/LAT data to 
investigate the empirical properties of the IC component of the SEDs of the FSRQs detected 
by Fermi. All the SED show apparent curvature and have a peak somewhere in the 10 MeV- 
1 GeV band. There is no correlation between the IC peak luminosity and energy for the 
sample of FSRQs detected by Fermi. The existence of such correlatio n has been claimed in 



the past for the lu minosity and the energy of the synchrotron peak (iGhisellini et al.lll998 
Fossati et al.lll998l ) for a sample of blazars (i.e. FSRQs and BL Lacs). Thus it might be that 



this correlation (if real) exists only when the two families of blazars are joined together and 
that any correlation for the FSRQs class is washed away by the presence of the additional 
EC component. Also the lack of correlation of the IC peak luminosity and frequency reveals 
that FSRQs are, unlike BL Lacs, part of a population with homogenous properties. 

We built average redshift-corrected SEDs in four different luminosity bins. The average 
SEDs are surprisingly similar as a function of luminosity (and redshift) as Fig. testifies. 
Approximating the SED with a power law with an index 2.4, while not producing the correct 
shape, allows the reader to compute a k-correction useful up to redshift ?a2. Beyond that 
this approximation is not valid. 



8.4. The Contribution to the Diffuse Background 



This work has important consequences for our understanding of the origins of the diffuse 



background. As pointed by several authors (e.g. Ilnoue &: Totanil 120091 ) and determined in 
this work, the spectrum of the diffuse emission arising from FSRQs shows curvature, due 
to the curved SEDs of these objects. We couple our model SED to our LF to predict the 
FSRQ contribution to the 10 keV to 100 GeV diffuse background. FSRQs produce ~10 % of 
the cosmic diffuse emission in the 1 MeV-10 GeV band. 
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Because of its good sensitivity Fermi has already resolved as much as 50 % of the total 
flux from FSRQs in the 100 MeV-10 GeV band. Our a nalysis indicates that the contribution 



of unresolved FSRQs to the IGRB flAbdo et aUbOlOdf ) is 9.66 



thE x io- 



ph cm 



sr 



and thus only 9.3^};q% (±3% systematic uncertainty) of the intensity of th e IGRB. This 
analysis is in good agreement with the results reported by lAbdo et al.l (J2010eJ) except above 
10 GeV where the use of a simple power law for the spectra of FSRQs was inadequate. 



Our results appear in reasonably good agreement with those of llnoue et al.l (120101 ) and 
of llnoue fc Totanil ( 20111 ). both in terms of spectral shape of the diffuse emission arising from 



FSRQs and its intensity. In their work, these authors rely on the sample of FSRQs and BL 
Lacs detected by EGRET. It is thus not surprising that their estimates of the contribution 
to the IGRB are slightly large r than ours. Fi nally, our estimate reported above is in good 
agreement with the results of iDermerl (120071 ) that predicted that FSRQs would produce 
~10-15 % of the 7-ray background. 

LAT-detected FSRQs represent only ~ 0.2 % of the parent population (see §[7j) and thus 
it is reasonable to expect that misaligned jets, although less luminous, but more numerous 
give a non-negligible contribution to the diffuse background. Our beaming model allowed us 
to explore this scenario. It is found that misaligned relativistic jets contribute ~30 % of the 
diffuse flux from the FSRQs class at redshift zero. If the Lorentz factor distribution does 
not change with redshift then the contribution of unresolved FSRQs and their misaligned 
siblings might be arou nd ~2.0xl0~ 6 ph cm -2 s" 1 and thus ~20% of the IGRB. Recently, 
Inoue &: Totanil (120 111) predicted that radio galaxies of both the FR-I and FR-II type might 
be able to account for ~25% of the intensity of the IGRB. In our work we found that 
FR-II alone could in principle (see above caveat) produce ~10% of the IGRB. It can be 
envisaged that once also the contribution to the IGRB of BL Lacs and their parents will be 
established, the total 7-ray emission from relativistic jets might account for some ~40-50 % 
of the intensity of the IGRB. 
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Table 4. Parameters of the beaming models described in the text. Parameters without an 
error estimate were kept fixed during the fitting stage. 



Paremeter 


Value 


Value 


k 


-2.03±0.70 


-2.43±0.11 


hi 


5.1±0.5 a 


5.0±0.5 b 


B 


3.04±0.08 


3.00±0.08 


7i 


5 


5 


72 


40 


40 


&i 


10 40 


10 38 




10 44 


10 42 


V 


4 


5 


X 2 /dof 


1.3 


1.5 



a In units of 10 23 . 
b In units of 10~ 26 . 
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A. Systematic Uncertainties 



The sources of systematic uncertainties in this analysis are: incompleteness (i.e. missing 
redshifts), detection efficiency, blazar variability, a nd EB L. A detailed discussion of some of 
these problems was already given in lAbdo et al.l feOlOd ) . Incompleteness in our sample is 
very small and introduces no appreciable systematic uncertainty as shown in § 14.31 



A.l. Detection Efficiency 



The detection efficiency used to deter mine the sky area surveyed by Fermi at any given 
flux is very important in this analysis fsee lAbdo et al.ll2010d. for a det ailed discussion). The 
detection efficiency used in this work was derived in lAbdo et al.l feOlOef ) under the assumption 
that the blazars spectra can be approximated by a power law. While this might be true over 
a small energy band, it becomes a problematic assumption over the full 100 MeV-100 GeV 
band covered by the LAT. In this Section we estimate directly the systematic uncertainties 
connected to thi s hypothesis. We performed 3 end-to-end Monte Carlo simulations of the 
Fermi sky (see lAbdo et al.ll2010d . for details), assigning randomly a curved spectrum to 
each source. These spectra are extracted from a library created using the ~100 observed 
spectra derived in § [5] varying the parameters of the measured spectra within their errors. 
The simulations were then analyzed to derive the detection efficiency reported in Fig. [161 In 
particular in order to detect a source a maximum likelihood fit with a power law spectrum is 
performed. This is done in order to reproduce the i nherent systematic uncertainty of fitting 
the curved spectrum of a source with a power law (lAbdo et al.ll2010ah . 



Fig. [16] shows the detection efficiency for a sample like that used in this analysis. Two 
aspects are noteworthy. First the efficiency at F 10 o = 10 -8 ph cm~ 2 s _1 (i.e. the lowest flux 
of this analysis by construction) is ~ 0.02 with a typical uncertainty of ±5 x 10~ 3 dictated 
by the small statistic in our simulations at the lowest fluxes. Second, at fluxes around 
Fioo ~ 10~ 7 ph cm' 2 s' 1 the detection efficiency becomes larger than 1.0. This effect is due 
to the fact that fitting a curved spectrum source with a power law yields t o an overestimate 
of the source flux by a factor ~ 10% (see also Fig. 8 in lAbdo et al.ll2010el ). Since Fig. [TBI is 
built as the ratio (in a given bin) of the number of sources detected with a given flux to the 
number of simulated sources with that flux, the effect mentioned above leads to a detection 
efficiency >1.0. 

In order to test the level of systematic uncertainty we derived the LF using the detec- 
tion efficiency reported in Fig. [TBI Given the "small" number of sources detected in the 3 
simulations, it was not possible to derive a two-dimensional detection efficiency as a func- 
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tion of flux and sp ectral index (like that used in § H] and derived for power law sources in 
Abdo et al.ll2010el ). For this reason the parameters of the distribution of photon indices of 
the FSRQ class cannot be derived from the analysis of the LF. As it is apparent from Table [3] 
most parameters of the LF derived in this section and those derived in § 14.21 are compatible 
within their statistical errors. The only parameter for which the difference is slightly larger 
than the statistical errors is a. The parameter a governs the trend of the redshift peak with 
luminosity and while its statistical error is in both case 0.03, the systematic error appears to 
be 0.05. This has very little impact on the analysis and the results of the previous sections 
are fully confirmed and robust against variations of the detection efficiency curve. As a 
further proof, the points of the de-evolved LF in Fig. [5] and Fig. [6] were computed using the 
detection efficiency of Fig. [16] while the shaded error region was computed using the model 
LF derived in § 14.21 that uses the detection efficiency for power law source. 

We performed an additional test, by shifting the detection efficiency curve in Fig. [16] to 
fluxes 10 % brighter than measured. The rightward shift is most dramatic as it increases the 
magnitude of the correction at faint fluxes. The shift is performed in order to account for 
uncertainties in the determination of the detection efficiency. The parameters of the LF are 
all consistent within statistical uncertainty with those found in this and the previous sections 
and reported in Table [3] The index of the low-luminosity slope of the LF becomes slightly 
steeper (i.e. 7i=0.47±0.18), and this yields a slightly larger contribution to the IGRB from 
FSRQs. We thus consider the typical systematic uncertainty connected to the estimate of 
the contribution to the IGRB to be ~3 % of the IGRB 100 MeV-100 GeV intensity. 



A. 2. Variability 



It is well known that blazars are inherently variable objects with variability in flux of 
up to a factor 10 or more. Throughout this work only average quantities (i.e. mean flux, 
mean luminosity and mean photon index) are used. 

It is not straightforward to determine how blazar variability affects the analysis pre- 
sented here. While the variabili ty pa tterns and amplitudes of blazars as a class are still not 
known both I Abdo et al.l ( ]2010d ) and I Abdo A., et al.J (120111 ) presented a detailed analysis of 
the variability of the brightest Fermi blazars. They report that the variability amplitude 
of the FSRQ class is generally larger than that of the BL Lac population. However, most 
sources (either bright or faint ones) exceed their average flux for less than 5-20 % of the 
monitored time (i.e. respectively 11 months or 2 years). This drastically reduces the possi- 
bility that FSRQs (or blazars more in general) are detected because of a few bright flaring 
episodes. The effect of high-amplitude variability connected with a rising density of sources 
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Fig. 16. — Detection efficiency as a function of flux for a population of sources with curved 
spectra similar to those of FSRQs determined in § 



at smaller flux might contaminate samples, as the one used here, with objects that formally 
should not be included. However, be cause of what reported above and the fl atness of the 
FSRQs source count distribution (see lAbdo et al.ll2010el ; lAbdo A., et al.Jl201lL and Fig. [2]) 
this effect is very likely marginal. 

Another, smaller, problem is connected to the dependence of the effective area on the 
direction of the incoming photon and the LAT detector fram^]. Short intense flares detected 
during favorable conditions (i.e. on axis and at azimuthal angles of ~0, ~90, ~180, or ~270 
degrees) might lead to a higher TS, increasing the likelihood of source detection. However, 
because of Fermi's con tinuous scanning o f the sky and because most flares are observed to 
last 10 days or longer (lAbdo et al.ll2010d ). the effect above has negligible influence on the 
analysis presented here. 

Finally, we believe that variability does not hamper the results of this analysis. Using 
even longer integration times (e.g. 2, or 3 years) will be the most efficient way to confirm 



7 Se e e.g. 



http://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm 



and 



Atwood et al. 
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the results of this analysis and dilute the effect of blazar variability. 

A. 3. Extragalactic Background Light 

Uncertainty in the level of the EBL, in particular at medium to high redshift, might 
in principle affect our analysis. Energetic 7-rays from FSRQs at high redshifts might be 
absorbed by the EBL and if this effect is not taken into account the source-frame luminosity 
would be underestimated. This would lead to wrong estimates of the space densities of 
FSRQs. However we believe this uncertainty is negligible. 

The uncertainty in the level of the EBL would impact the estimate of the k-correction 
which allows us to determine the source-frame luminosities. As shown in Fig. [IOJ neglecting 
the EBL at once and adopting a simple power-law spectrum with a photon index of 2.4 
(instead of the average SED determined in § \5§ introduces an uncertainty of <4 % on the 
value of the k-correction at z=3. Since all the Fermi FSRQs are detected within this redshift, 
this uncertainty produces no appreciable impact on the determination of the luminosity 
function. 
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